/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package simulatepushpull;
import java.util.Random;
/**
 *
 * @author vipR
 */
public class Gamma {
    private Random RNG;
    public Gamma(Random RNG)
    {
        this.RNG = RNG;
    }
     public double sampleGamma(double k, double theta) 
     {
    boolean accept = false;
    if (k < 1) {
 // Weibull algorithm
 double c = (1 / k);
 double d = ((1 - k) * Math.pow(k, (k / (1 - k))));
 double u, v, z, e, x;
 do {
  u = RNG.nextDouble();
  v = RNG.nextDouble();
  z = -Math.log(u);
  e = -Math.log(v);
  x = Math.pow(z, c);
  if ((z + e) >= (d + x)) {
   accept = true;
  }
 } while (!accept);
 return (x * theta);
    } else {
 // Cheng's algorithm
 double b = (k - Math.log(4));
 double c = (k + Math.sqrt(2 * k - 1));
 double lam = Math.sqrt(2 * k - 1);
 double cheng = (1 + Math.log(4.5));
 double u, v, x, y, z, r;
 do {
  u = RNG.nextDouble();
  v = RNG.nextDouble();
  y = ((1 / lam) * Math.log(v / (1 - v)));
  x = (k * Math.exp(y));
  z = (u * v * v);
  r = (b + (c * y) - x);
  if ((r >= ((4.5 * z) - cheng)) ||
                    (r >= Math.log(z))) {
   accept = true;
  }
 } while (!accept);
 return (x * theta);
    }
  }

    
}
